Role of dipolar correlations in the IR spectra of water and ice 
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Abstract 

We report simulated infrared spectra of deuterated water and ice using Car-Parrinello molecu- 
lar dynamics with maximally localized Wannier functions. Experimental features are accurately 
reproduced within the harmonic approximation. By decomposing the lineshapes in terms of intra 
and intermolecular dipole correlation functions we find that short-range intermolecular dynamic 
charge fluctuations associated to hydrogen bonds are prominent over the entire spectral range. Our 
analysis reveals the origin of several spectral features and identifies network bending modes in the 
far IR range. 

PACS numbers: 36.20.Ng, 71.15.Pd, 32.10.Dk, 78.20.Ci 
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I. INTRODUCTION 



It is well known that correlations among the dipole moments of hydrogen bonded 
molecules greatly enhance the static dielectric response of water and ice (see e.g. Ref. [l| 
and references therein). However, the role played by dipolar correlations in the dynamic 
dielectric response of hydrogen bonded systems is not well understood, although it has often 
been speculated that intermolecular couplings affect significantly infrared (IR) and Raman 
spectra . A difficulty in modeling these spectra is that the dynamic dielectric sus- 

ceptibility depends on the adiabatic response of the electrons to nuclear dynamics. Another 
difficulty, in principle even greater, is that nuclear quantum effects cannot be discarded a 
priori in water systems at, or below, room temperature, and it is beyond current theoretical 
and computational capabilities to fully account for dynamic quantum effects in systems of 
this complexity. Recent studies have shown, however, that the IR absorption lineshapes of 
liquid water and ice are overall reproduced remarkably well by ab initio molecular dynamics 
(AIMD) simulations, which treat nuclear dynamics classically but include explicitly, and 
accurately, the quantum adiabatic response of the electrons 7, f|. 

In this paper we adopt the AIMD framework and use maximally localized Wannier func- 
tions (MLWF) to separate intra and inter molecular contributions to the dynamic dipolar 
correlations in liquid water and ice. We find that hydrogen bonds (H-bonds) make short- 
range intermolecular fluctuations of the electronic charge as important as the intramolecular 
fluctuations over the entire spectral range. Our analysis sheds new light into the origin of 
several spectral features and provides an unambiguous identification of network bending 



modes in the far IR region [9|. 

The IR absorption coefficient per unit length a(uj) is related to the refractive index n(uj) 
and the imaginary part of the dielectric constant e"(uS) by a{uj)n{uj) = (lo / c)e" (uj) . Within 
linear response theory, a(uj) is g iven by the power spectrum of the time correlation function 
of the total dipole operator [lOj. Following common practice, we approximate the quantum 
time correlation function with the classical one, i.e. with (M(0)M(t)), where M is the total 
dipole moment in the simulation cell and the brackets indicate classical ensemble average. 
Since there is only one classical correlation function while the quantum time correlation 
function can be expressed in several equivalent ways, this substitution leads to formulae 
for the IR absorption coefficient characterized by different prefactors known as quantum 
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correction factors ll). We adopt here the so-called harmonic approximation (HA) which 
follows by replacing the Kubo-transformed quantum correlation function with the classical 
one, leading to: 

a(u)n(u) = I dte- wt {M{0)M{t)) (1) 



3cV 



where (3 = jknT) 1 is the inverse temperature. In the harmonic regime HA is exact 12]. 



Ramirez et al. 



ll| showed that HA is the only correction factor that satisfies the fluctuation- 



dissipation theorem in addition to detailed balance. The same authors also found that HA 
performs better than the other quantum correction factors for one dimensional anharmonic 
potentials that model different H-bond scenarios. 

II. COMPUTATIONAL DETAILS 



In this work, we compute the IR spectra of D2O ice and water using a Car-Parrinello 
(CP) scheme 13], in which maximally localized Wannier functions (MLWFs) [yj] represent 
"on the fly" the electronic wave functions 15j. In the simulation of water (ice), we use a 
periodically repeated cubic (orthorhombic) cell with 64 (96) D2O molecules at the density 
of 1.1 (1.0) g/ml. For ice we use a proton disordered ice Ih configuration generated as 
in Ref. jl]. We adopt norm-conserving pseudopotentials r1 an d the DFT Perdew-Burke- 
Ernzerhof (PBE) functional for exch ang e and correlation [171 ]. using a plane- wave cutoff of 
60 Ry for water and of 85 Ry for ice [18j. We use an integration time step of 7 a.u. (0.17 fs) 
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131 ] . The temperature is 
with mass Qn = 1 X 10 6 . 



and a fictitious electron mass of 350 a.u. in the CP equations 
controlled by coupling the nuclei to a single Nose thermostat [] 
In this way the average temperature of the system is set to 330 K in water and to 268 K in 
ice. 



At each time step we assign a dipole Hi = r l D + r % D2 + 6r l — 2 J2i=i,4 r Wi to the i-th. 
molecule in the cell. Here r l x are the positions of the nuclei (X = Di, D2, O) and of the 
four MLWF centers (X = Wi, with I = 1,4) associated to the eight valence electrons of a 



water molecule 



20| . The i-th dipole is conventionally located at the molecular center of 



mass r\ and depends on the local environment. The total dipole moment is M = J2i=i n fa 
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221 ] and the IR spectra are computed with Eq. (CQ). The time correlation function is 



averaged over an equilibrated MD trajectory of 10 ps in ice and of 23 ps in water. A 
Gaussian window function 1 231 is used in the discrete Fourier transform. 
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III. RESULTS AND DISCUSSION 



A. Calculated IR spectra 



The calculated spectra are compared with experiment in Fig. [T] 24] . As in previous HA 
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FIG. 1: (color online). Experimental (upper panels) and calculated IR s 



Experimental data are from Ref. 25[ for H2O ice at 100K, and from Ref. 



aectra (lower panels). 



2d] for D 2 ice at 150K. 



D2O water data are from Ref. |27( at 295K and from Ref. [9j at 293K. Experimental data for D2O 
ice are not available over the entire range. 

calculations [?], 0] the overall agreement between theory and experiment is good. For ice, 
the positions of the band maxima corresponding to H-bond stretching, libration, bending, 
combination, and OD stretching are: 200 (222), 660 (640), 1140 (1210), 1660 (1650), and 
2120 (2425) cm -1 respectively (wavenumbers in parentheses are experimental values from 



Ref. 



28j). For water, the corresponding values are: 200 (187), 510 (505), 1140 (1215), 1550 



(1555), and 2160 (2450) cm -1 respectively. The details of the experimental features are well 
reproduced, e.g. the asymmetric shape of the OD stretching bands, the skewed ice libration 
band, and even the small combination bands. The good performance of the harmonic ap- 
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proximation may appear somewhat surprising given the large shifts and broadenings of the 
vibrational frequencies of an isolated molecule that arise in condensed phase. Furthermore, 
significant anharmonicity has been detected in the excited hydrogen stretching modes in liq- 



uid water 



291 ] . IR absorption, however, probes equilibrium properties which are dominated 



by the vibrational ground-state in the hydrogen stretching region. Then, the overall similar- 
ity between calculated and observed spectra suggests that the main effects of anharmonicity 
are sufficiently well captured by classical dynamics. 

The largest discrepancy between simulation and experiment occurs in the hydrogen 
stretching modes, which are redshifted compared to experiment by ~300 cm -1 . Part of this 
error originates from our choice of the mass \x for the fictitious dynamics of the electrons. We 
can quantify this effect by calculating the IR stretching band in ice from Born-Oppenheimer 
(BO) molecular dynamics simulations, in which the electrons are kept in the instantaneous 
ground-state by minimizing the energy functional at each time step. In CP simulations 
this condition is enforced dynamically and the outcome depends on the choice of \x. BO 



calculations show that a redshi 



t of up to ~80 cm in the IR stretching band of ice can be 



attributed to our choice of \i 30f] . However, even when the BO separation is strictly enforced 
the IR stretching band in ice is more than 200 cm -1 below experiment. The adopted DFT 
approximation is a likely cause of this error. For instance, in the H2O water monomer, PBE 
yields stretching frequencies that are ~140 cm -1 lower than the corresponding experimental 
values 32j. Taking the isotopic mass effect into account we should expect a redshift of ~100 
cm -1 for the stretching band of D2O in gas phase. The actual effect that we observe in 
condensed phase, ~200 cm -1 , is larger than that and is consistent with the known H-bond 
over binding present in PBE water. Other sources of inaccuracy in the calculated spectra 
include the finite basis set, the cell size and the effect of temperature. Quantum anharmonic 
corrections are beyond our approach. 

The similarity between the spectra of solid and liquid water is consistent with the standard 
picture of a tetrahedral, ice-like local structure of the liquid. With more disorder in the 
liquid phase, modes close in frequency overlap with each other, leading to broader bands 
than in the solid. Fluctuating local fields and coupling between neighboring molecules cause 
additional broadening. The relative importance of local fields and intermolecular couplings 
can be quantified in terms of intra and inter spectral contributions. 
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B. Intra vs. inter contribution 

While separating intra and intermolecular contributions would be difficult or even im- 
possible in experiment, in theory it can be easily achieved due to the additive nature of 
two-body correlation functions. In Fig. [2] we report the intra and inter contributions to 
the IR spectra of ice and water. The total spectrum is computed from Eq. (pQ), in which 
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FIG. 2: (color online). Intra (red dash) and inter (blue dot) contributions to the calculated IR 
spectra (black solid). Lower panels magnify the (0, 800) cm^ 1 region. The inset shows the power 
spectrum of the autocorrelation function of the molecular angle #(DOD), and of the H-bond angle 
#(0...0...0), in arbitrary units. The arrows indicate positions in cm . 

(M(0)M(t)) = (Y,ij fJ>i(0)iJ,j(t)) . The intra contribution corresponds to retaining only the 
terms with i = j, the remaining terms add up to the inter contribution. 

In water, the most prominent band around 2100 cm -1 is due to the overlap of the sym- 
metric and asymmetric OD stretching modes together with the overtone of the bending 
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mode 



281 ]. With two maxima, the inter contribution has larger intensity and lower fre- 



quency compared to its intra counterpart. This gives rise to the three broad features and the 
asymmetric shape of the band observed in experiments. Similar effects occur in ice, where 



three main features have been experimentally identified in the stretching band [4j,[26j. In this 
case, however, we can not identify distinct intermolecular features within our chosen Fourier 

n 

filter |24J]. Two peaks resembling the experimental features appear in the inter contribution 
when we reduce the filtering but they are too closely separated compared to experiment. In 
contrast to the stretching modes, the inter contributions to the bending modes around 1140 
cm -1 are negative in both ice and water. This effect has a simple intuitive explanation: due 
to the electrostatic attraction between opposite charges, OD stretching modes are correlated 
between neighboring molecules whereas bending modes are anti-correlated. In addition, we 
find that the broad shape of the small combination band at ~1550 cm -1 in water (~1660 
cm -1 in ice) is mostly due to the inter contribution. 

The far infrared region is more complex. Two features have been identified in the libration 
band of the Raman spectrum of the liquid, but only one mode was believed to be IR active 
on the basis of a simple tetrahedral model with C2 V symmetry 28J. A more recent IR 
experiment, however, suggested that both modes are visible in the liquid These modes 
have their counterpart in ice and, indeed, we find two IR features at ~460 cm -1 and at 
~660 cm -1 in our calculated spectrum. Both features have intra and inter contributions 
but, interestingly, the inter contribution at ~460 cm -1 is negative, largely cancelling its 
intra counterpart. Instead the inter contribution at 660 cm -1 is small and positive, slightly 
enhancing the corresponding spectral feature. As a result, the ice libration band has the 
skewed shape observed in experiments (Fig. [1]). In the liquid, the situation is similar with 
the features shifted to ~370 cm -1 and ~590 cm -1 . Broadening due to disorder makes the 
lineshape look more symmetric. 

The features at ~50 cm -1 in water and at ~200 cm -1 in water and ice are H-bond 
network modes with bending and stretching character, respectively In Ref. Q, the 

inter character of the 200 cm -1 mode was emphasized. A more accurate analysis reveals, 
however, that both intra and inter contributions are present, with an inter/intra ratio of 2:1 
in ice and of 1:1 in water. This analysis is facilitated in the case of ice, where the network 
modes are well separated from the libration band. 



7 



C. Identification of the network bending mode 

Interestingly, we also understand the origin of the feature at ~50 cm -1 . Albeit weak, this 
feature has been identified in experiments particularly below room temperature [3|. In our 
simulation it appears as a weak intramolecular feature in both water and ice, where the peak 
is at ~70 cm" 1 . Its intra character is surprising because the feature originates from network 
bending modes, as shown by the power spectrum of the autocorrelation function of the angle 
between two adjacent H-bonds, i.e. the angle defined by the oxygen atoms of three molecules 
linked by H-bonds. The corresponding spectrum for ice is reported in the inset of Fig. [2] and 
is sharply peaked at ~70 cm -1 . We also report in the same inset the power spectrum of the 
autocorrelation function of the angle between the two covalent OD bonds in a molecule. As 
expected this spectrum is sharply peaked at ~1140 cm" 1 but, interestingly, it also exhibits 
a weaker feature at ~70 cm -1 due to the modulation of the molecular angle induced by the 
network bending modes. Thus, the absence of inter character in the IR spectrum derives 
from the negligible intermolecular charge fluctuation associated to network bending modes, 
while the intra contribution reflects the modulation of the molecular dipole moment induced 
by these modes. 



D. Spatial extent of dynamic dipolar correlations 

To gain insight on the spatial extent of the intermolecular correlations, we extend the 
approach of Ref. [l| to the dynamic domain. We define the density of the intermolecular 
dipole correlation function as: 

D inter (r,t) = — 5a^< E E/*(°Vi(t)fy(*)> (2) 

i=l,N j^i 

where r is the distance between the center of mass of molecule j at time t and molecule i at 
time 0. Sij(t) = 1 if at time t molecule j is in spherical shell (r, r + Ar) centered on molecule 
i at time 0, whereas = otherwise. The corresponding integrated contribution to the 
IR spectrum up to distance R is given by: 

a(u)n(u)\^%=—-f dr dte~ mt D inter {r,t) (3) 

OCV JO J~oo 

The R dependence of several peak intensities is reported in Fig. [3J The curves indicate 
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FIG. 3: Inter contributions of the peak intensities vs. R (Eq. ([3])). The maximum value of R 
corresponds to half the size of the simulation cell. The points at 15 bohr report the peak values 
from Fig. (2J which have been calculated without spherical cutoff. 



that the most important intermolecular correlations for all the modes occur when R is 
comprised between 4 and 6 bohr, i.e. in correspondence with the first coordination shell, 
in both water and ice. Beyond 6 bohr only minor additional contributions are present, 
similar to the findings for static correlations in Ref. Short-range dynamic correlations 
are a direct manifestation of the presence of H-bonds between adjacent molecules. While 
static dipolar correlations due to H-bonds always enhance the dielectric response, dynamic 
correlations do enhance the system response at some frequencies but they also suppress it 
at other frequencies. 



IV. CONCLUSION 



Given the importance of intermolecular correlations, only models that include their effect 
at a sufficient level of detail can realistically describe the observed experimental features. 
An important result of our study is that these correlations play an equally important role 
over the entire spectral range and are not limited to the network modes that characterize 
the far IR range. We expect that similar dynamic correlations should also be crucial to 
interpret other spectroscopic data, such as Raman, which probe the adiabatic dynamics of 
the electrons. It would be difficult to include properly these effects in simulations based 
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on classical empirical potentials. Our analysis shows that liquid water and ice are systems 
in which the molecules are strongly correlated, suggesting that their IR spectra cannot be 
simply explained in terms of a weakly perturbed collection of individual molecules. Finally, 
the decomposition technique introduced in this paper is not restricted to bulk water and ice, 
but can also be applied to other systems and to water at interfaces. In the latter case one 
may expect that characteristic signatures of hydrophilic and hydrophobic interfaces should 
emerge from this analysis. 
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